1. Field of the Invention
The present invention relates to a method and apparatus for analyzing chemical systems. More particularly, the present invention relates to a method and apparatus for analyzing large molecules using approximate molecular geometry.
2. Description of the Related Art
In classical mechanics, it is assumed that particles can possess arbitrary amounts of energy. However, on a molecular scale, this assumption fails. Therefore, quantum mechanics was developed to predict the amount of energy possessed by a particle and the probability that a particle will be positioned at a certain point. Quantum mechanics are explained in detail in texts such as Physical Chemistry, which is hereby incorporated by reference in its entirety for all purposes.
Quantum mechanics can be used to study chemical compounds and ions which are molecules which have acquired a charge due to the gain or loss of electron(s). Hereinafter, the terms "chemical system" and "chemical compound" are used synonymously to encompass chemical compounds and ions.
As part of the quantum mechanical analysis, it is necessary to solve Schrodinger's equation which is represented as: EQU H.psi.=E.psi.
where H is the Hamiltonian operator which is a grouping of parameters for the particular chemical system; .psi. is the wave function (an operator) which replaces the classical concept of trajectory; and E is the energy of the particle.
According to the Born approximation, the probability density of a particle (the probability that a particle will occupy a certain position) is given by the square of the wavefunction, .psi..sup.2.
In order to solve Schrodinger's equation it is necessary to determine the electron-electron interaction for each of the electrons of the molecule (in each molecular orbital). The potential energy of each electron depends upon the potential energy of every other electron. Thus, the potential energies are initially guessed and solving Schrodinger's equation is preformed iteratively. The cyclic calculations are performed until the orbitals and the energies obtained are insignificantly different from those used at the start of the latest cycle. The solutions are then said to be self-consistent, and accepted as solutions of the problem.
In a many electron system, an analytical solution is impossible. However, very powerful computational techniques are available to provide detailed and reliable numerical solutions for the wavefunctions and the energies.
A solution to Schrodinger's equation allows one to optimize molecular structure by minimizing the energy of the particles. There are at least two computational techniques to determine molecular structure including the semi-empirical methods, which combine empirical data and quantum mechanical calculations, and ab initio methods, which involve calculation of orbital integrals based on first principles.
An example of a semi-empirical method is AMI, PM3, and MNDO which are found in MOPAC 93 (Molecular Orbital PACkage). To understand the MOPAC 93 method, it is helpful to look at an example. The exact calculations of this example are described on pages 118-123 of the MOPAC 93 Manual, which is hereby incorporated by reference in its entirety for all purposes. The steps which are described below are shown in the flow chart of FIG. 1.
The system to be examined is a regular hexagon of hydrogen atoms in which the H--H distance is 0.98316. Of course, a regular hexagon of hydrogen atoms is not a stable system; the only reason it is being used is to demonstrate the working of a self consistent field ("SCF") calculation.
In the analysis, first an approximate atomic distance is input (step 1 in FIG. 1) which is converted to an interatomic distance matrix, shown below, where the distance is given in Angstroms.
______________________________________ Interatomic Distance Matrix (.ANG.) Atom 1 2 3 4 5 6 ______________________________________ 1 0.0000 2 0.9832 0.0000 3 1.7029 0.9832 0.0000 4 1.9663 1.7029 0.9832 0.0000 5 1.7029 1.9663 1.7029 0.9832 0.0000 6 0.9832 1.7029 1.9663 1.7029 0.9832 0.0000 ______________________________________
Then, the interatomic distance matrix is used to determine a one-electron matrix (step 2 in FIG. 1) which ignores the effects of other electrons. The one-electron matrix is shown below:
__________________________________________________________________________ One-electron matrix (eV) Atom 1 2 3 4 5 6 __________________________________________________________________________ 1 -51.7124 2 -3.2457 -51.7124 3 -1.0970 -3.2457 -51.7124 4 -0.6992 -1.0970 -3.2457 -51.7124 5 -1.0970 -0.6992 -1.0970 -3.2457 -51.7124 6 -3.2457 -1.0970 -0.6992 -1.0970 -3.2457 -51.7124 __________________________________________________________________________
The analog of the one-electron matrix in Huckel theory is the Alphas and Betas. See section 16.4(a) of P. W. Atkins, Physical Chemistry, Third Edition (1986) for a discussion on Huckel's theory. The diagonals are the energy an electron would have if it were only on one atom (the wavefunction of one atom).
The off-diagonals are the energy the electron would have if it were in the wavefunction of two different atoms. As can be seen, as the distance between the atoms becomes greater (refer to the distance matrix), the energy of the electron becomes smaller (is inversely proportional to the distance). This is because there is a larger interaction between close atoms.
Then, inter-electron interactions are considered. A calculation is made as to the energy a first electron would have if it were in the field of a second electron (step 3 in FIG. 1). The resulting two electron integrals are shown below, where the energy is given in electron volts:
______________________________________ Two-Electron Integrals (eV) Atom 1 2 3 4 5 6 ______________________________________ 1 12.8480 2 9.6585 12.8480 3 7.0635 9.6585 12.8480 4 6.3622 7.0732 9.6585 12.8480 5 7.0635 6.3622 7.0732 9.6585 12.8480 6 9.6585 7.0635 6.3622 7.0732 9.6585 12.8480 ______________________________________
The energy is positive because the electrons are repelling one another.
A density matrix is necessary in order to calculate the Fock matrix, but, in turn, the Fock matrix is necessary in order to calculate the density matrix. To break this impasse, a guessed density matrix is formed (step 4 in FIG. 1). The guess is very crude: all off-diagonal matrix elements are set to zero, and all on diagonal terms on any atom are set equal to the core charge of that atom divided by the number of atomic orbitals. In other words, all electrons are assumed to be localized on their respective atoms. Our starting guess for H.sub.6 consists of a unit matrix (1's on the diagonals). The starting density matrix is not shown.
Each iteration of the SCF calculation consists of assembling a Fock matrix from the one-electron matrix, the two-electron integrals, and the density matrix, diagonalizing it to obtain the eigenvectors, and finally reassembling the density matrix. At some point the change in the density matrix drops below a preset limit. When this happens we say that the field is self-consistent. These steps will now be carried out for the H.sub.6 system.
Formation of the initial Fock matrix (step 5 in FIG. 1) is simple, as there are no off-diagonal terms in the density matrix. Therefore the off-diagonal terms for the initial Fock matrix are the same as those terms in the one-electron matrix.
The on-diagonal terms in the initial Fock matrix Faa are modified by the electrostatic field of all the electrons in the system except the electron or fraction of an electron in the atomic orbital .phi..sub.a. Consider F(1,1). The total initial population of .phi..sub.1 is 1.0 (from the starting density matrix), composed of equal amounts of .alpha. and .beta. electron density. An electron in .phi..sub.1 would therefore experience the electrostatic repulsion of half an electron. An electron cannot repel itself; however, it will be repelled by its partner electron of opposite spin.
In addition, each electron will be affected, normally repelled, by the electrostatic field of all the electrons on all the other atoms. Each atom has one electron, so the total energy of an electron, i.e., the diagonal initial Fock matrix element, is given by: EQU F(1,1)=-51.7124+1/2(12.848)+2(9.6585+7.0635)+6.3622.
The initial Fock matrix is obtained by adding the two-electron terms to the one-electron matrix. The elements of the initial Fock matrix represent the sum of the one and two electron interactions. For the system of six hydrogen atoms, this has the following form:
__________________________________________________________________________ Initial Fock Matrix (eV) Atom 1 2 3 4 5 6 __________________________________________________________________________ 1 -5.4823 2 -3.2457 -5.4823 3 -1.0970 -3.2457 -5.4823 4 -0.6992 -1.0970 -3.2457 -5.4823 5 -1.0970 -0.6992 -1.0970 -3.2457 -5.4823 6 -3.2457 -1.0970 -6.6992 -1.0970 -3.2457 -5.4823 __________________________________________________________________________
The initial Fock matrix is then diagonalized (step 6 in FIG. 1) to yield the following set of eigenvalues, or one-electron energies, and eigenvectors, or molecular orbitals. The energy levels listed to the left represent the molecular orbitals. Because the example is directed to an H.sub.6 system, only 6 electrons are present and only the lowest three energy levels are occupied. To the left of the energy levels are the molecular orbital ("MO") coefficients.
__________________________________________________________________________ Molecular Orbital Coefficients Energy Level 1 2 3 4 5 6 __________________________________________________________________________ 6 -0.4857 0.4082 -0.4082 0.4082 -0.4082 0.4082 -0.4082 5 -1.8388 0.5774 -0.2887 -0.2887 0.5774 -0.2887 -0.2887 4 -1.8388 0.0000 0.5000 -0.5000 0.0000 0.5000 -0.5000 3 -6.9317 0.5774 0.2887 -0.2887 -0.5774 -0.2887 0.2887 2 -6.9317 0.0000 0.5000 0.5000 0.0000 -0.5000 -0.5000 1 -14.8670 0.4082 0.4082 0.4082 0.4082 0.4082 0.4082 __________________________________________________________________________
A new density matrix, shown below, is then calculated (step 7 in FIG. 1) from the diagonalized Fock matrix:
______________________________________ Density Matrix (eV) Atom 1 2 3 4 5 6 ______________________________________ 1 1.0000 2 0.6667 1.0000 3 0.0000 0.6667 1.0000 4 -0.3333 0.0000 0.6667 1.0000 5 0.0000 -0.3333 0.0000 0.6667 1.0000 6 0.6667 0.0000 -0.3333 0.0000 0.6667 1.0000 ______________________________________
The diagonal terms of the new density matrix are the sum of the squares of the MO coefficients for the first three energy levels times two. For example, the 1--1 one diagonal of the new density matrix is the MO 1 coefficient for energy level 1, 0.4082, squared, plus the MO 1 coefficient for energy level 2, 0.0000, squared, plus the MO 1 coefficient for energy level 3, 0.5774, squared. The sum of these three terms is then multiplied by two, because the spins of the electrons are in opposite directions.
The off-diagonal terms are determined in the same manner except the two different MO coefficients are multiplied, rather than squaring one MO coefficient. For example, the atom 1-2 term is determined by multiplying MO coefficients 1 and 2 in energy level one, multiplying MO coefficients 1 and 2 for energy level 2 and multiplying MO coefficients 1 and 2 for energy level three. Then these three products are added together and the sum is doubled.
A new Fock matrix can then be constructed (step 8 in FIG. 1) using the new density matrix. Each element of the new Fock matrix (below) is calculated by taking a corresponding term from the one-electron matrix and subtracting from it one half the product of the corresponding density matrix term and the corresponding two electron integral term.
__________________________________________________________________________ Second Fock Matrix (eV) Atom 1 2 3 4 5 6 __________________________________________________________________________ 1 -5.4823 2 -6.4652 -5.4823 3 -1.0970 -6.4652 -5.4823 4 +0.3611 -1.0970 -6.4652 -5.4823 5 -1.0970 +0.3611 -1.0970 -6.4652 -5.4823 6 -6.4652 -1.0970 +0.3611 -1.0970 -6.4652 -5.4823 __________________________________________________________________________
Because the on-diagonal terms of the new density matrix are one, the on-diagonal terms of the new Fock matrix are identical to those in the initial Fock matrix. The off-diagonal terms are however changed. The off-diagonal terms are modified to allow for exchange interactions. (Note that not all exchange terms are stabilizing.)
The new Fock matrix element F(1,2), for example, is evaluated as follows: EQU f(1,2)=-3.2457-1/2(0.6667) (9.6583) eV
Next, it is determined if there is self consistency (step 9 in FIG. 1). This is done by diagonalization of the new Fock matrix and comparing the total electron energy of the resulting matrix with the matrix originally obtained. Alternatively, the new density matrix could be used to determine if there is self consistency. Other elements of the analysis could also be used. In this case the same set of eigenvectors are obtained as were initially obtained (the eigenvalues obviously will be different). Thus, because there is no change, the comparison has shown that the difference is below a preset threshold and there is self-consistency, and the SCF calculation is complete (10 in FIG. 1). This is quite unusual and the simplicity is one reason the H.sub.6 hexagonal system was chosen. In general, several iterations are necessary in order to obtain an SCF. If there is not self consistency, the Fock matrix is diagonalized (step 6 in FIG. 1) and the steps subsequent thereto are repeated.
Currently quantum chemical analysis on large chemical systems such as enzymes is impractical even with the fastest semi-empirical methods such as MOPAC 93 and ab initio methods, such as GAUSSIAN, ATMOL, GRADSET, GAMES and HONDO. However, MOPAC 93 and all other semi-empirical and ab initio methods involve diagonalization of the Fock matrix. The computation time of diagonalization rises by the cube of the size of the matrix (number of atoms in the molecule). Even on the fastest supercomputer large molecules, such as enzymes, would require many years to process. Mathematicians have tried for many years to quicken the diagonalization step calculations. However, the analysis time for the diagonalization step still rises as the third power as the size of the Fock matrix.